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NASA Technical Memorandum X -62000 
DIGITAL FILTER SYNTHESIS PROGRAM 
By Robert M. Munoz and Richard A. Moyer 

Page 6 - Equation 20 should read: 

0 (Z) Bi Z-i 0 (Z) B 2 Z-2 I (Z) A 0 I (Z) Ai Z _1 I (Z) A 2 Z“ 2 

o (z) + ^ + + 

JDq -£>0 -^O -DO -DO 

Page 6 - Equation 21 should read: 

0 (Z) = 0 (Z) [-Bi Z" 1 -B 2 Z" 2 ] + I (Z) [Ao + A[ Z" 1 + A 2 Z" 2 ] 

Page 12 - example problem load deck second card from the end should read: 

l. 12.5663706 10. .001 

Page 13 - Appendix B 

Example problem print -out third line should read: 

1.0 12.57 10.00 .00100 
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TECHNICAL MEMORANDUM X-62000 


DIGITAL FILTER SYNTHESIS PROGRAM 

By Robert M. Munoz and Richard A. Moyer 

Ames Research Center 
Moffett Field, Calif. 9^035 

ABSTRACT 


In order to perform computations on the differential equations of con- 
tinuous time varying quantities using the general purpose digital computer, it 
is necessary to form the difference equations or computational algorithms that 
approximate the differential equations. To do this, the "Digital Filter Syn- 
thesis Program" has-been written. This program written in Fortran IV uses the 
bilinear transform method to approximate linear differential equations with 
constant coefficients . It allows inputs to be presented as functions of "s" 
(complex frequency) in either the factored or unfactored form and it presents 
the coefficients of the difference equations as intermediate outputs . The 
program tabulates input and output data lists representing the discrete form 
of time varying input and output quantities for the transfer functions repre- 
sented by the differential equations. A subroutine for graphical display of 
the outputs by use of the print plot technique is also included in the program. 
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DIGITAL FILTER SYNTHESIS PROGRAM 
By Robert M. Munoz and Richard A. Moyer 


SUMMARY 

In this report, a discussion of the theory and operation of a digital 
filter synthesis program is given. This program was written specifically to 
aid in the computer simulation studies of space instrument systems but it has 
broad general application. It allows any continuous function of a complex 
variable to be expressed in approximate form as a computational algorithm or 
difference equation . 


INTRODUCTION 


Digital filtering has received a great amount of attention in recent 
years and a number of methods have evolved for digital filter synthesis. 

Rader and Gold (ref. l) have given a good summary of many of these methods and 
Kuo and Kaiser (ref. 2) have dedicated a chapter of their book to the s ame 
topics. In this report, a brief description will be given of the automatic 
implementation of one method of digital filter synthesis originated by 
Steiglitz (ref. 3) and discussed by Rader and Gold. This bilinear transform 
method allows synthesis of the difference equation for digital filtering with 
the initial specifications for filter performance given in terms of the 
analog prototype in the frequency domain. In other words, the digital filter 
can be derived from F(s) the transfer function for the real time analog 
equivalent network. The pertinent mathematical results are also summarized 
briefly. 

The digital filter synthesis program can perform the operation of digital 
filtering on any input data once the difference equation has been developed. 
Since the theory underlying synthesis of the analog prototype is well devel- 
oped (see refs. 4 and 5 ), one can design a digital filter with any desired 
frequency or transient response by the use of this program. 

The functions performed by the digital filter synthesis program are 
illustrated in figure 1. The first step, however, which is not a part of this 
program, is the synthesis of the analog prototype of the filter. The proto- 
type is a filter configuration scaled for a frequency of 1 radian per second 
and 1 ohm impedance level. As an example of a prototype filter configuration, 
figure 2 shows the transfer function and the frequency response function for a 
Butterworth filter, a filter that has a second order maximally flat amplitude 
response. Any prototype filter form is possible; however, care should be 
exercised in using filters of fifth order and above because numerical insta- 
bility may result from the cumulative effect of round off errors in such 



filters. Where higher order filters are required , they should he produced hy 
a cascade of lower order filters developed from the factored form of the 
prototype transfer function. High pass, low pass, and hand pass 
configurations as well as combinations are acceptable. 

The digital filter synthesis program is written in Fortran IV and oper- 
ates under the IBM 7040/709^ DCS (Direct Couple System) with the IBSYS 
executive monitor. Plotting is available for a maximum of 500 points per 
variable. The number of points plotted by the program may be any fraction 
less than the number of points computed. This fraction is determined by the 
plot ratio and the total number of points equals run time divided by the 
input sampling interval. 


BILINEAR TRANSFORM METHOD 


The bilinear transform method -is an approximate method because digital 
data only approximately represent analog data and because the bilinear trans- 
form warps the frequency scale . Where very many samples of a given analog 
waveform are taken in an interval large with respect to filter time constants, 
which is the mode of operation expected, a very good approximation is 
obtained. However, the program is useful even when this is not the case, but 
care should be taken to compute or calibrate the performance of the filter 
when fewer than 10 samples are used in the interval corresponding to the 
shortest filter time constant. Methods for compensating the frequency warping 
effect are available but have not been incorporated in the program. 

A bilinear transformation is the process of forming all the coefficients 
of the digital difference equation which controls the actual digital filter- 
ing. The form of the difference equation is such that the current data out- 
put sample is formed from the present data input, past data input, and past 
data output samples only. Other types of digital filter operations using 
future data input samples are not possible in this program. The structures 
of filters of this type are discussed at length by Arabadjis (ref. 6) and 
others . 


DESIGN PROCEDURE 


The design procedure used in this program starts with the choice of a 

sampling interval t. This is the time in seconds between each entry in the 

digital input data. Once this time interval is known and the digital filter 

critical frequencies have been established, it is possible to compute the 

equivalent analog filter critical frequencies according to the following 
formula: 


OJa .T 

“ai “ tan ~*t 


( 1 ) 


2 
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where the are the fundamental quantities used in frequency scaling for 

the analog prototype, w&i are the corresponding quantities for the digital 
filter, T is the interval between data points, and i is an index. 

Given a prototype filter expressed in the following form: 

ya s n 

Filter function = — — (2) 

2b m s 

where a and b are arbitrary constants, m and n are real positive integers, 
and s is the complex frequency variable, the program computes a new 
frequency-scaled function as follows: 


F(s) 


• n 

La s 
n 




m 


(3) 


where the primes represent scaled quantities . The next step in the synthesis 
process is the transformation from the s to the z plane by means of the 
bilinear transformation approximation. 


Choice of Sampling Interval 

The sampling frequency should be large in comparison to the frequencies 
of interest and especially the critical frequencies for reasons mentioned 
earlier. In order to simplify the discussion, let us continue using the 
example of figure 2 and carry out the numerical aspects of the computations 
as we go along. The data we will be dealing with in this example are in the 
order of cycles per second, we will choose 1000 samples of data per second. 
This is more than necessary but should make the approximation to the analog 
prototype quite good; namely, 


T = 0.001 




Choice of the 3 dB Cutoff Frequency as a Critical Frequency 

To pursue the numerical example further, we assign the 3 dB point for the 
Butterworth filter to be a critical frequency at 2 Hz: 

= 2?!^ = 12.56 rad/sec (5) 

where f^ is the desired 3 dB cutoff frequency in Hz for the digital filter. 
Then u a is computed as follows: 

w a = - « 0.00628 ( 6) 
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Choice of a Filter Form 


The prototype of figure 2 has "been chosen for this example and the analog 
transfer function is repeated here: 

F(s) = i (7) 

s^ + \[ 2 . s + 1 


Calculation of the Frequency-Scaled Parameters for the Prototype 

For synthesis on the real frequency axis we may substitute jw for s 
in equation (7) and, wherever ju> occurs, scale this variable by the fre- 
quency scaling ratio which, for this case, is equal to w a as follows: 



F*(u>) = i 

(M) 

(8) 

F(s) can then be 

calculated by substituting back as follows: 



F* ( s ) = 

S 2 s/28 n 

— _ + + 1 

w a 2 w a 

(9) 


F»( s ) = i 

b' 2 s 2 + b’ is + b' 0 

(10) 

Coefficients are 

matched to give the b values where 



b 1 2 = -i- - 2.5^XL0 4 
w a a 

(11) 


b»i = — = 2 . 25XL0 2 

(12) 

*’o = 1 

and can be computed for any order filter by similar procedures. Thus, 

(13) 


F'(s) = — 

(2.5^X10 4 )s 2 + (2.25X10 2 )s +1 

(1*0 
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Bilinear Transformation of F(s) 
How to calculate F(z) replace s "by z - l/z + 1 

F( z ) = i 


(2-54X10 4 ) - i + (2.25X10 2 ) - v + 1 

v z+1 ' z+1 


z - 1 


z+1 

Multiplying numerator and denominator by (z + l)^ gives: 

z 2 + 2z +1 


F(z) = 


2.54KL0 4 z 2 - 5-08XL0 4 z + 2.54X10 4 + 2.25X10 2 z 2 - 2.25X10 2 + 1 


2 2 


( 15 ) 


(16) 


Multiplying numerator and denominator again by z~ 2 gives 


F(z) 


1 + 2z 1 + z 2 


25625 - 508OOZ- 1 + 25176z- 2 


(17) 


Since F(z) represents the digital equivalent of a transfer function, it can 
be represented by the following expression: 


F( z ) = - AnZ ^ = °M 


SB n z -rl ^ 


(18) 


where 

0(z) output data function of z 
l(z) input data function of z 

Ah numerator coefficients of the digital transfer function 

Bn denominator coefficients of the digital transfer function 

The function 0(z) represents the desired output digital data list and we can 
solve for this quantity by simply rearranging equation (l8) as follows: 


0(z) = 


I(z)A q + l(z)AiZ 1 + l(z)A 2 z 


-2 


B q + B 1 Z -1 + B 2 z 


7-2. 


(19) 


Dividing through by B 0 and simplifying yields 
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0(z) = - 


l(z)A2Z -2 


( 20 ) 


0(z)BiZ -1 0(z)B 2 z 2 1(z)A q + l(z)AiZ -1 + 

B 0 B 0 B 0 B 0 B 0 

Since each term on the right-hand side of equation (20) is divided by B Q , a 
new set of constant coefficients is defined for simplicity as follows: 

0(z) = 0(z)[-B'xZ - 1 - B' 2 z - 2] + l(z) [A ’ 0 + A'jZ - 1 + A'z - 2] (21) 


This is the difference equation that represents the nearest digital equivalent 
to the analog prototype within the limitations allowed by the input sampling 
intervals for the example. The following is a list of these coefficients: 

A * 0 = 3 . 913 x 10“ 5 1 


A* i = 7 . 826x10" 5 
a' 2 - 3.913XLO- 5 > 
B'l = 1.982 
B* 2 = O.98238 


( 22 ) 


The program lists these coefficients for each filter configuration. 


GRAPHICAL OUTPUTS 


One of the special features of this program is the graphical displays of 
data lists that are used as inputs and outputs of the data lists that are 
used as inputs and outputs of the filters. Graphs are obtained in the print 
out by the method of plotting that was developed at the University of 
Michigan. Special characters are used to form a reference grid and then 
variable points from the same listing machine are used to generate tabulated 
digital data. The graphs formed by this method are only accurate to approxi- 
mately ±1 percent, but they do allow an almost instantaneous graphical output 
to be generated, thereby providing improved man-machine communications. Other 
methods such as off-line plotting have the advantage of greater accuracy but 
the disadvantage of long turn-around time. In trouble-shooting and debugging 
operations, this turn-around time is cumulative and under some circumstances 
can be devastating. Some of the limitations in accuracy of the list plotting 
method used in this program are compensated for by simultaneously listing the 
digital data so that if it is desired to measure a point more accurately than 
is possible on the graph, the tabulated data will give the result to any 
degree of precision required. It is not necessary to plot each point in 
input or output data. A plot ratio allows one out of every n points to be 
plotted where n is a positive integer. For example, if the plot ratio is 
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10 , computations will be performed and the required accuracy will be achieved 
with all input digital data points though only 1 out of every 10 points will 
be plotted. 


SUMMARY OF SUBROUTINES 


In this program, a number of subroutines that accomplish the functions 
indicated in figure 1 are discussed here. The BILIN routine obtains filtering 
constants from a transfer function. ISCALE is a routine that sets up the time 
scale abscissa and number of samples for the graphing routines. The FILTER 
subroutine takes one input and gives one output for each entry, retaining only 
those past values which it needs. PLT subroutine will scale the ordinate 
values for the plot and give the first two calls to the UMPLOT subroutine 
(PL0T1 and PL0T2). The R array of the program is a collection of all ordinate 
points for the graph with LR = the number of such points. Any combination or 
arrangement of graphs of the different operations and filters is possible. 

The plot is then completed by successive PL0T3 calls (.1-5.) for each operation 
and a PLOT^ call. A plot exceeding 500 points is made possible only by 
changing the image dimension of subroutine PLT. The following is a listing 
of the subroutines with their arguments and definitions . 


COEF (N,C) 
ST1002 


XDEN (N,S,B) 
ST1003 


Stores binomial coefficients. The coefficients of the x 2 , 

x 3 , . . . , J terms are placed in C(3) , C(4) , . . . , 

C(N + l) . C(l) and C(2) are fixed by the programmer. 

Assumes coefficient of s^ to be 1. The coefficients of 
sN-1, s® -2 , . . . , s constant are assumed to be in S(l) , 
S(2) , . . , , S(N - l) , S(N) , respectively. XDEN makes the 
substitution (z - l) N , (z - l) N+1 (z + l) , (z - l) N_2 (z+l) 2 , 
. . . , (z - l)(z + l)^ - -'- for s^ _ l, s-^“2, . . . , s and 
multiplies the constant by (z + l) . The coefficients of 
z j , j = 0, 1, 2, . . . , N are given in B(l) , B(2.) , • • •, 
B(N), B(N + l) • This subroutine is part of the bilinear 
transformation . 


FILTER (VAL,YY, 
A,B,X,Y,NN,M,NNR) 
ST1004 
VAL 
YY 

A*s and B’s 
X's and Y's 

M 

NN 

NNR 


Performs filter function 


Xi - input data list l(z) 

Yi - output data list 0(z) 

Filtering coefficients supplied by BILIN 

Arrays containing past Input and output quantities needed 
for the calculation of YY. 

Order of the filter 
Equal to M + 1 
Number of the sample 
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PSMPY (P,NP,Q, 
NQ,RR,NRR) 

ST1005 

PROD (EK,A,B,MU, 
M,P,NP) 

STIOO 6 

EK 

A 

B 

MJ 

M 

P 

UP 


RCOEF (H,SS,S) 
ST1007 
N 
S 

XINP (VAL,.TI,K) 
ST1008 

TI 

K 

NOTE: 


PLT (KL,HE,R,IB, 
KLL,IFL) 

ST1009 

KL 

RT 

R 

LR 

KLL 

IFL 


Multiplies polynomials P (an array with, dimension NP) and 
Q (an array with dimension NQ) giving the resultant poly- 
nomial in array RR with dimension NRR 

Produces coefficients of polynomials from its roots. 


Scaling factor (the coefficient of the highest term in 
degree; 1. for this program) 

Real part of roots (an array) 

Imaginary part of roots (an array) 

Array specifying number of times each root is repeated (all 
set to 1.) 

Number of different roots 

Array of coefficients starting with constant term rearranged 
in subroutine RCOEF. 

Degree of resulting polynomial 

A pair of conjugate roots is counted as one root . The 
conjugate is assumed where the root location has an 
imaginary component . 

Generates coefficients of polynomial S from complex roots 
in SS 

The number of roots 

Real coefficients of resulting polynomial 

Generates filter inputs and is a subroutine supplied by the 
programmer. A single input is calculated on each entry and 
placed in VAL . 

The sampling time interval 
The nth input point 

The programmer should compile this subroutine for each run 
to insure proper inputs to the filter. 

Provides the first two calls (PL0T1 and PL0T2) to the plot 
routine UMPLOT. 

All arguments to this subroutine must be calculated and 
supplied to it. 

Number of samples to be plotted 

An array of times (abscissa) for the plot 

An array containing all the ordinate values for scaling 

purposes . 

Number of elements in array R 
Number of abscissa grid lines 

Flag; positive formal entry; negative to call PL0T2 only. 
This subroutine scales the ordinate values with the help of 
subroutines RANGE and SCALE. SCALE is a library subroutine. 
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NOTE: 


ISCALE (TIME,TI , 
IW,NP,KLL) 

ST1011 

NP 

KLL 

time 

TI 

IW 

RANGE (RMIN, 
RMINI, LS) 

STIOIO 

RMIN 

RMINI 

LS 


SCALE (RNG , 

rmin,io.,sf) 

RNG 

RMIN 

10 . 

SF 

RNG 

SCALE 


BILIN (LD,LN,FR, 
TI,A,B,C0N,IFL) 
ST1012 
LD 

LN 

FR 

TI 

A and B 


The abscissa and ordinate are interchanged so that the 
abscissa may be several pages long. This necessitates the 
change in sign for time (RT) values and this is accomplished 
by the subroutine automatically. 


Calculates arguments NP and KLL from the other three 
arguments . 


Total number of points in the sample 
Number of abscissa grid lines 
Length of run in seconds 
Sampling period in seconds 

Plot ratio, total number points/number points plotted 
Calculates a rounded minimum for scaling. 


minimum value to be plotted, floating point 

Minimum value to be plotted, integer 

Number of shifts 

Enter with least value in RMIN 

Exit with new rounded value in RMIN 

Examples: 


Enter RMIN 


Exit RMIN 


0 .014+ 
0.0025+ • 
1.4 
-0 .013 


0.0100— 
0 . 0020 — 
1.0 
- 0.020 


Calculates a maximum value 


Maximum value to be plotted 
Minimum 

Available inches minus 1 
scale factor 

ll.xSF + RMIN, calculates new value of maximum based on 
11 inches . 

is a NASA Ames utility subroutine written by J. A. Jeske 
available at this time on the disk. No deck card is needed. 

Calculates coefficients for the difference equation 


Degree of the polynomial or number of poles in the 
denominator of the transfer function f(s) 

Degree of the polynomial or number of roots in the numerator 
of the transfer function „f(s) 

The frequency scaling ratio 
Sampling period in seconds 

Arrays containing coefficients for the difference equation 
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C0N 

IEL 

Positive 

Negative 


Leading constant of the transfer function 
Input flag 

f(s) function in terms of coefficients of polynomials 
f(s) function in terms of poles and zeros of f(s) 

All arguments are inputs to the subroutine except arrays 
A and B. See flow chart and inputs to BILIN 


The following plot subroutines are NASA Ames modifications from SHAKE 
library routine UMPLOT. 


PL0T1 (NSC ALE, 

NHL ,NSBH ,NBL ,NSBV ) 


NSCALE 
NSCALE(l) = 

NSCALE (2) = 
NSCALE(3) = 
NSCALE(4) = 

NSCALE(5) = 

NHL 

NSBH 

NBL 

NSBV 


An array of dimension 5 that supplies the subroutine with 

grid and scale factor information 

0, standard grid and scale factors (see note (a)) 

0, grid and scale factors are as defined in NSCALE(2)- 
NSCALE(5) 

1, scale factor such that printed values of the ordinate 
are 10^ times the actual values 

J , j digits will appear to the right of the decimal point 
in printed ordinate values (j < 8) 

K, scale factor such that printed values of the abscissa 
are 10^ times the actual values 

M, M digits will appear to the right of the decimal point 
in printed abscissa values (M < 8) 

The number of horizontal grid lines (NHL > 0) 

The number of spaces between horizontal grid lines 
(NSBH > 0) 

The number of vertical grid lines (NBL > 0) 

The number of spaces between vertical grid lines 
(NSBV > 0, and NSBV*NBL < 119 ) 


NOTE (a): 


Standard scale factors correspond to values of I, J, K, 
and M of 0, 3, 0, 3, respectively. 


PL0T2 (IMAGE, 

XMAX,XMIN,YMAX, 

XMIN,IDIM) 

IMAGE 

XMAX 

XMIN 

XMAX 

YMIN 

IDIM 


An array, dimensioned IDIM, •which is used as the image 
region for the plot 

The value of the abscissa at the right most grid line 

The value of the abscissa at the left most grid line 

The value of the ordinate at the upper most grid line 

The value of the ordinate at the lower most grid line 

The dimension of the array IMAGE, where IDIM is at least 
equal to N*(NSBH*NHL + l) whera^. , pi-w/pl'l 

n = [ (k/6 ) + (1/2) + (i/2)(-ir K 1_2LK/2J) ] 

and where K = NSBV*NBL + 1 

(The square brackets in the formula for N signify 
"integral value.”) 
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PL0T3 (1HX,RX,RT,L) 


IHX 

X is the character used for plotting the points whose 
abscissas are stored in array RX and whose ordinates 
are stored in array RT 

RX 

The array containing the abscissas of the points to be 
plotted . 

RT 

The array containing the ordinates of the points to be 
plotted 

L 

The number of points to be plotted 

PL0T4 (20,20Hbbbb 

This subroutine initiates the plot and provides the 

Tbbb IbbbMbbbEbbb ) 

ordinate label 


FLOW CHARTS 


Flow charts for the main program and for the subroutine BILIN are given 
in figures 3 and 4. 


EXAMPLE 


A load deck for the example filter given in figure 2 with a sinusoidal 
input is shown in appendix A. Note that the input function which may be any 
analytic or data input function is entered in subroutine XINP to be compiled 
each time the problem is run. A list of input values X, output values Y, 
and the graphical outputs of the program is also given for this example for 
a time varying between 0 and 10 seconds. The A’s and B's listed are the 
coefficients of the difference equation for the digital filter. 


INPUT FORMAT 


Tables I and II describe the input format for the Digital Filter Syn- 
thesis Program. This format contains the information necessary to produce the 
coefficients of the difference equation as well as ancillary inputs to 
determine the character of the input data lists, output data lists, run t ime , 
and plot characteristics . 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 9^035, Sept. 25, 1967 
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EXAMPLE PROBLEM LOAD DECK 
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LPT ( M )~ 
LPB(M)- 



LPT(M) =LPT( M ) *KK 
LPB(M)=LPB(M)*LK 
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T = (T)0 
X=(I )d 



P { 2 ) =K N 
Q( 2 ) = K 

IF(K.EQ.N) GO TO 
CALL COEF (KN,P) 
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IF(NR.LT.NN) RETURN 
DO 40 1 = 1 , NRX 
XU )=X (1 + 1 ) 

Y ( I ) =Y ( 1 + 1 ) 



40 CONTINUE 
LAST - 1013 
RETURN 
END 
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C PROD LIST »REF»DECK»M94»XR7 

CPROD (F4) T417 W.M.SYN (2/18/64) 

C... CALCULATE THE COEFFICIENTS OF A POLYNOMIAL FROM ITS ROOTS 



DIMENSION A(l)> B(l), MU(1)» P(l)» R(3) 
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END 

$IBFTC ST 1007 

SUBROUTINE RCOEF ( N f SS, S) 

DIMENSION FIRST! 1) 

DIMENSION S ( 50 ) t D ( 50 ) ,DI ( 50 ) » MU ( 50 ) t B ( 50 ) 



COMPLEX SSI 50) 

DATA M0/50*l/ 

DO 10 1=1, N 
D ( I ) = REAL ( SS ( I ) ) 
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RNG = 0.0 
DO 130 1=1, LR 
IF (R < I ) .GT.RNG ) RNG=R( I ) 
IF(R( I) .LT.RMIN) RMIN=R t I ) 



130 CONTINUE 

CALL RANGE (RMIN»RMINI*LS) 
CALL SCALE(RNG,RMIN,10.,SF) 
RNG=11.*SF+RMIN 
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SUBROUTINE I S G A L E ( T I M E » T I » I W , N P , K L L ) 

NP IS THE NUMBER OF POINTS IN THE EXPERIMENT 
NP-TIME/T I 
ID=IW*5 

KLL IS NUMBER OF GRID LINES FOR ABSCISSA 



KLL - NP/ID 
IC = ID*KLL 
IF(NP.GT.IC) GO TO 50 
RETURN 
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I F ( K ) 1 1 » 1 2 1 9 
CONTINUE 
C ( 1 ) = 1. 

C ( 2 ) = FLOAT ( K) 



CALL COEF (K,C ) 

CALL PSMPY(C,Kf AtLNvClt IORD) 
NK - IORD + 1 
DO 10 1=1, NK 
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TABLE I— INPUTS FOR DIGITAL FILTER SYNTHESIS PROGRAM 


Comment card to be printed verbatim (columns 1-80) 

Column 

1-5 LD Degree of the polynomial in the denom- 

inator of the filter transfer function 
f(s). Exception is noted with con- 
jugate roots (poles of F(s)). 

6-10 LN Degree of the polynomial in the numer- 

ator of the transformation function 
f(s). Exception is noted with conju- 
gate roots (zeros of F(s)). LN is 
zero if the numerator is a constant. 

11-15 IW Plot ratio; ratio of the total number 

of samples to number of plotted samples 

16-20 JJ Flag for print/plot options 

Print and plot 
Plot only 
Print only 
Do neither 

21-25 IFL Input flag 

positive f(s) function in terms of coefficients 
of polynomials 

negative ' f(s) function in terms of poles and 
zeros of f(s) . 


3 

1-10 

CON 

Leading constant of the transfer 
function f(s) 

(all floating 
point) 

11-20 

FR 

The scaling frequency ratio 


21-30 

TIME 

Total time of the filter run in seconds 


31-40 

TI 

Time interval between successive points 


in the input data in seconds 

The next card (or cards if required) provide inputs to subroutine BHIN in 
either the factored or unfactored form of the prototype transfer function 
depending on whether IFL is positive or negative. The fields are 10 columns 
of 8 floating point numbers per card, as illustrated in table II. The 
denominator of the transfer function is given first and the numerator follows 
immediately. The coefficients of these polynomials must be entered as real 



Card 

1 

2 

ill integers 
ght justified) 
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numbers . Roots of the polynomials entered in factored form must be complex 
numbers. If a root has a nonzero imaginary part, its conjugate is assumed to 
be another root and must not be entered as input, that is, complex conjugate 
pairs are always assumed and one complex number represents the pair. 
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TABLE II.- SAMPLE INPUTS FOR SUBROUTINE BILIN 


1 , 


F(s) 


s 2 + 3s + 8 
s 3 + 3s 2 + 3s + 1 


LN = 2 
LD = 3 
IFL = 1 
C0N - 1 


Floating point mode 

Column 1 11 21 30 

3- 3- 1* 

3 - 8 . 


3 LN = 0 

2. F(s) = - LD = 1 

(s + 1.5 - J2.4)(s + 1.5 - j2.4) IFL = -1 

C0N = 3 


Complex mode 

Column 1 11 20 

- 1.5 2.4 


F(s ) = Ms + 1)(2 - 1 + 33)(s - 1 - J3) 
(s + 2)(s + 3 + Jl)(s + 3 - jl) 


Complex mode 
Column 1 11 21 

-2. 0.0 -3. 

— 1 . 0.0 1 . 


31 to 

1 . 

3 - 


4 . 


l.5(s 2 + 2s + 3) 

F(s) = 

s 10 +3s 9 +8s 8 +7s 7 -+2s 6 +2s 5 +s 4 +3s 3 +4s 2 +s+6 


LN = 2 
LD = 2 
IFL = -1 
C0N = 4 


LN = 2 
LD = 10 
IFL = 1 
C0N =1.5 


Column 1 11 

3. 8. 

1 . 6 . 

2- 3- 


Floating point mode 

21 ' 31 to 51 

7. 2. 2. 1. 


61 

3- 


71 80 

4. 
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f AGS BLANK NOT FILMED. 


Digital 

input 

data 



- Functional block diagram of the digital filter synthesis program 


Figure 1 . 





Figure 3.- Flow chart 


main program 







TRANSFORM NUMERATOR | 
ONLY, CALL XDEN ONCE 
NUMERATOR — (Z+l) LD I 
CALL CO.E.F 


TRANSFORM NUMERATOR AND 
DENOMINATOR CALL XDEN TWICE 

s LD , s LDH ._-s°~-(z-i) | - D ,{z-o LO ''(z+r)-.^-{z+n LD 

| t LN >T LN-I_^_ T o_^ (Z „ 0 LN >(Z _ 0 LN-I( Z+ ,) (Z+1)LN 

K IS DIFFERENCE IN DEGREE OF THE POLYNOMIALS 
NUMERATOR AND DENOMINATOR OF F(S) 


DEN — DEN*(Z+I)-K 
CALL CO.EF 


CALCULATE 
ARRAYS A SB 


NUM — «-NUM*| 
(Z+I) K 
CALL COEF 


C RETURN ) 


Figure 4.- Flow chart for subroutine BILIN . 
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